Using Fano factors to determine certain types of gene autoregulation

The expression of one gene might be regulated by its corresponding protein, which is called autoregulation. Although gene regulation is a central topic in biology, autoregulation is much less studied. In general, it is extremely difficult to determine the existence of autoregulation with direct biochemical approaches. Nevertheless, some papers have observed that certain types of autoregulations are linked to noise levels in gene expression. We generalize these results by two propositions on discrete-state continuous-time Markov chains. These two propositions form a simple but robust method to infer the existence of autoregulation in certain scenarios from gene expression data. This method only depends on the Fano factor, namely the ratio of variance and mean of the gene expression level. Compared to other methods for inferring autoregulation, our method only requires non-interventional one-time data, and does not need to estimate parameters. Besides, our method has few restrictions on the model. We apply this method to four groups of experimental data and find some genes that might have autoregulation. Some inferred autoregulations have been verified by experiments or other theoretical works.


PKC
: : Figure 1: An example of a GRN in human T cells [91]. Each vertex is a gene. Each arrow is a regulatory relationship. Notice that it has no directed cycle.

Introduction
In general, genes are transcribed to mRNAs and then translated to proteins. We can use the abundance of mRNA or protein to represent the expression levels of genes. Both the synthesis and degradation of mRNAs/proteins can be affected (activated or inhibited) by the expression levels of other genes [42], which is called (mutual) gene regulation. Genes and their regulatory relations form a gene regulatory network (GRN) [16], generally represented as a directed graph: each vertex is a gene, and each directed edge is a regulatory relationship. See Fig. 1 for an example of a GRN. The expression of one gene could promote/repress its own expression, which is called positive/negative autoregulation [11]. Autoregulation is very common in E. coli [63]. Positive autoregulation is also called autocatalysis or autoactivation, and negative autoregulation is also called autorepression [4,21]. For instance, HOX proteins form and maintain spatially inhomogeneous expression of HOX genes [64]. For genes with position-specific expressions during development, it is common that the increase of one gene can further increase or decrease its level [81]. Autoregulation has the effect of stabilizing transposons in genomes [6], which can affect cell behavior [41,79]. Autoregulation can also stabilize the cell phenotype [2], which is related to cancer development [96,52,15,90].
While countless works infer the regulatory relationships between different genes (the GRN structure) [86], determining the existence of autoregulation is an equally important yet less-studied field. Due to technical limitations, it is difficult and sometimes impossible to directly detect autoregulation in experiments. Instead, we can measure gene expression profiles and infer the existence of autoregulation. In this paper, we consider a specific data type: measure the expression levels of certain genes without intervention for a single cell (which reaches stationarity) at a single time point, and repeat for many different cells to obtain a probability distribution for expression levels. Such single-cell non-interventional one-time gene expression data can be obtained with a relatively low cost [48].
With such single-cell level data for one gene V , we can calculate the ratio of variance and mean of the expression level (mRNA or protein count). This quantity is called the variance-to-mean ratio (VMR) or the Fano factor. Many papers that study gene expression systems with autoregulations have found that negative autoregulation can decrease noise (smaller VMR), and positive autoregulation can increase noise (larger VMR) [70,68,30,51,25,18,17]. This means VMR can be used to infer the existence of autoregulation. We generalize the above observation and develop two mathematical results that use VMR to determine the existence of autoregulation. They apply to some genes that have autoregulation. For genes without autoregulation, these results cannot determine that autoregulation does not exist. We apply these results to four experimental gene expression data sets and detect some genes that might have autoregulation.
We start with some setup and introduce our main results (Section 2). Then we cite some previous works on this topic and compare them with our results (Section 3). For a single gene that is not regulated by other genes (Section 4) and multiple genes that regulate each other (Section 5), we develop mathematical results to identify the existence of autoregulation. These two mathematical sections can be skipped. We summarize the procedure of our method and apply it to experimental data (Section 6). We finish with some conclusions and discussions (Section 7).

Setup and main results
One possible mechanism of "the increase of one gene's expression level further increases its expression level" is a positive feedback loop between two genes [31]. Here V 1 and V 2 promote each other, so that the increase of V 1 increases V 2 , which in return further increases V 1 . We should not regard this feedback loop as autoregulation. When we define autoregulation for a gene V , we should fix environmental factors and other genes that regulate V , and observe whether the expression level of V can affect itself. If V is in a feedback loop that contains other genes, then those genes (which regulate V and are regulated by V ) cannot be fixed when we change V . Therefore, it is essentially difficult to determine whether V has autoregulation in this scenario. In the following, we need to assume that V is not contained in a feedback loop that involves other genes.
The actual gene expression mechanism might be complicated. Besides other genes/factors that can regulate a gene, for a gene V itself, it might switch between inactivated (off) and activated (on) states [10]. These states correspond to different transcription rates to produce mRNAs. When mRNAs are translated into proteins, those proteins might affect the transition of gene activation states, which forms autoregulation [23]. See Fig. 2 for an illustration. Therefore, for a gene V , we should regard the gene activation state, mRNA count, and protein count as a triplet of random variables (G, M, P ), which depend on each other.
When we fix environmental factors and other genes that affect V , the triplet (G, M, P ) should follow a continuous-time Markov chain. A possible state is the gene activation inactivated gene Figure 2: The mechanism of gene expression. A gene might switch between inactivated state and activated state, which correspond to different transcription rates. Gene is transcribed into mRNAs, which are translated into proteins. Proteins might (auto)regulate the state transition of the corresponding gene. state on/off (for G), the mRNA count on Z (for M ), and the protein count on Z (for P ). Thus the total space is {0, 1} × Z × Z. When we consider the expression level M or P (but have no access to the value of G), sometimes itself is Markovian (its dynamics can be fully determined by itself, without the knowledge of G), and we call this scenario "autonomous". In other cases, M or P itself is no longer Markovian (its dynamics explicitly depends on G), and we call this scenario "non-autonomous". We need to consider the triplet (G, M, P ) in the non-autonomous scenario. This is similar to a hidden Markov model, where a two-dimensional Markov chain is no longer Markovian if projected to one dimension (since this dimension depends on the other dimension).
For the autonomous scenario, we can fully classify autoregulation for a gene V . Assume environmental factors and other genes that affect the expression of V are kept at constants. Define the expression level (mRNA count for example) of one cell to be X = n, the mRNA synthesis rate at X = n − 1 to be f n , and the degradation rate for each mRNA molecule at X = n to be g n . This is a standard continuous-time Markov chain on Z with transition rates 1 ∆t P[X(t + ∆t) = n | X(t) = n − 1] = f n , Define the relative growth rate h n = f n /g n . If there is no autoregulation, then h n is a constant. Positive autoregulation means h n > h n−1 for some n, so that f n > f n−1 and/or g n < g n−1 ; negative autoregulation means h n < h n−1 for some n, so that f n < f n−1 and/or g n > g n−1 . Notice that we can have h n > h n−1 for some n and h n ′ < h n ′ −1 for some other n ′ , meaning that positive autoregulation and negative autoregulation can both exist for the same gene, but occur at different expression levels. Such non-monotonicity in regulating gene expression often appear in reality [1]. For the non-autonomous scenario, we can still define autoregulation. Consider the expression level X of V (mRNA count or protein count) and its interior factor I. If X is the mRNA count, then I is the gene state; if X is the protein count, then I is the gene state and the mRNA count. If there is no autoregulation, then X cannot affect I, and for each value of I, the relative growth rate h n of X is a constant. If X can affect I, or h n is not a constant, then there is autoregulation. When X can affect I, there is a directed cycle (X → I → X), and the change of X can affect itself through I. In this case, it is not always easy to distinguish between positive autoregulation and negative autoregulation.
Quantitatively, for the autonomous scenario, when we fix other factors that might regulate this gene V , if V has no autoregulation, then h n = f n /g n is a constant h for all n. In this case, the stationary distribution of V satisfies P(X = n)/P(X = n − 1) = h/n, meaning that the distribution is Poissonian with parameter h, P(X = n) = h n e −h /n!, and VMR = 1. If there exists positive autoregulation of certain forms, VMR > 1; if there exists negative autoregulation of certain forms, VMR < 1. However, such results are derived by assuming that f n , g n take certain functional forms, such as linear functions [54,57], quadratic functions [24], or Hill functions [67]. There are other papers that consider Markov chain models in gene expression/regulation [32,61,65,14,62,45], but the role of VMR is not thoroughly studied.
In this paper, we generalize the above result of inferring autoregulation with VMR by dropping the restrictions on parameters. Consider a gene V in a known GRN, and assume it is not regulated by other genes, or assume other factors that regulate V are fixed. Assume we have the autonomous scenario, meaning that its expression level X = n satisfies a general Markov chain with synthesis rate f n and per molecule degradation rate g n . We do not add any restrictions on f n and g n . Use the single-cell non-interventional one-time gene expression data to calculate the VMR of V . Proposition 1 states that VMR > 1 or VMR < 1 means the existence of positive/negative autoregulation.
Nevertheless, the autonomous condition requires some assumptions, and often does not hold in reality [5,39,33,32]. Consider a gene V that is not regulated by other genes, and has no autoregulation. The mRNA count or the protein count is regulated by the gene activation state (an interior factor), which cannot be fixed. Due to this non-controllable factor, there might be transcriptional bursting [60,19] or translational bursting [8], where transcription or translation can occur in bursts, and we have VMR > 1. This does not mean that Proposition 1 is wrong. Instead, it means that the expression level itself is not Markovian, and the scenario is non-autonomous. In this scenario, we should apply Proposition 2, described below, which states that no autoregulation means VMR ≥ 1.
We extend the idea of inferring autoregulation with VMR to a gene that is regulated by other genes, or with non-autonomous expression. Consider a gene V ′ in a known GRN. Denote other genes that regulate V ′ and the interior factors (gene state and/or mRNA count) of V ′ by F . Denote the values of V ′ , F as X, Y . Assume V ′ is not contained in a feedback loop, and assume g n , the per molecule degradation rate of V ′ at X = n, is not regulated by other genes or its interior factors (gene state and/or mRNA count). We do not add any restrictions on the synthesis rate f n . Proposition 2 states that if V ′ has no autoregulation, then VMR(X) ≥ 1. Therefore, VMR(X) < 1 means autoregulation for V ′ . Proposition 2 is derived in a "one-step" Markov chain model, where at one time point, only transitions to the nearest neighbors are allowed: (X = n, Y = a) → (X = n + 1, Y = a), (X = n, Y = a) → (X = n−1, Y = a), and (X = n, Y = a) → (X = n, Y = a ′ ). This one-step Markov chain model is the most common approach in stochastic representations of gene regulation [70,30,54,51,17]. Recently, there are some studies that consider "multistep" Markov chain models, where at one time point, the change of mRNA/protein count can be accompanied with the change of other factors, such as the gene state [7,43,74]. For example, the following transition is allowed: (G * , M = n) → (G, M = n + 1). In this multi-step model, Proposition 2 is no longer valid: even without autoregulation, it is possible that VMR(X) < 1. Consider an example that the production of one mRNA molecule needs many steps of gene state transition, and the gene returns to the initial step after producing one mRNA molecule: Since there are many steps, the total time for one cycle of G 1 → · · · → G k → G 1 + M can be highly deterministic, such as 1 second. Assume the degradation probability for each mRNA molecule in 1 second is 0.01. Then the mRNA count is highly concentrated near 100, and VMR(X) < 1 (close to 0.5 in numerical simulations).
Since multi-step models allow more transitions, they are more general than one-step models. However, it is still a question that whether such generalizations are necessary, since one-step models have good fitting with experimental data [38,18,9]. Proposition 2 provides a method to verify this problem: If a gene has VMR(X) < 1, but we use other methods to determine that it has no autoregulation, then Proposition 2 states that onestep models deviate from reality, and multi-step models should be adopted. Therefore, when one-step models hold, Proposition 2 is a valid method to determine the existence of autoregulation; when one-step models do not hold, combined with other methods to determine autoregulation, Proposition 2 can detect the failure of one-step models.
In the scenario that Proposition 2 may apply, if VMR ≥ 1, Proposition 2 cannot determine whether autoregulation exists. In fact, with VMR, or even the full probability distribution, we might not distinguish a non-autonomous system with autoregulation from a non-autonomous system without autoregulation, which both have VMR ≥ 1 [9]. In the non-autonomous scenario, we only focus on the less complicated case of VMR < 1, and derive Proposition 2 that firmly links VMR and autoregulation.
In reality, Proposition 1 and Proposition 2 can only apply to a few genes (which are not regulated by other genes or have VMR < 1), and they cannot determine negative results. Thus the inference results about autoregulation are a few "yes" and many "we do not know". Besides, for the results inferred by Proposition 1, especially those with VMR > 1 (positive autoregulation), we cannot verify whether their expression is autonomous, and the inference results are less reliable.
Current experimental methods can hardly determine the existence of autoregulation, and to determine that a gene does not have autoregulation is even more difficult. Therefore, about whether genes in a GRN have autoregulation, experimentally, we do not have "yes" or "no", but a few "yes" and many "we do not know". Thus there is no gold standard to thoroughly evaluate the performance of our inference results. We can only report that some genes inferred by our method to have autoregulation are also verified by experiments or other inference methods to have autoregulation. Besides, if the result by Proposition 2 does not match with other methods, it is possible that the one-step model fails. Instead, in Section A, we test our methods with numerical simulations, and the performances of both Propositions are satisfactory.

Related works
There have been some results of inferring autoregulation with VMR [70,68,30,51,25,18,17]. However, these VMR-based methods have various restrictions on the model, and some of them are derived by applying linear noise approximations, which are not always reliable in gene regulatory networks [72].
Besides VMR-based methods, there are other mathematical approaches to infer the existence of autoregulation in gene expression [59,92,22,73,37,97,35,34]. We introduce some works and compare them with our method. (A) Sanchez-Castillo et al. [59] considered an autoregressive model for multiple genes. This method (1) needs time series data; (2) requires the dynamics to be linear; (3) estimates a group of parameters. (B) Xing et al. [92] applied causal inference to a complicated gene expression model. This method (1) needs promoter sequences and information on transcription factor binding sites; (2) requires linearity for certain steps; (3) estimates a group of parameters. (C) Feigelman et al. [22] applied a Bayesian method for model selection. This method (1) needs time series data; (2) estimates a group of parameters. (D) Veerman et al. [73] considered the probability-generating function of a propagator model. This method (1) needs time series data; (2) estimates a group of parameters; (3) needs to approximate a Cauchy integral. (E) Jia et al. [37] compared the relaxation rate with degradation rate. This method (1) needs interventional data; (2) only works for a single gene that is not regulated by other genes; (3) requires that the per molecule degradation rate is a constant.
Compared to other more complicated methods, VMR-based methods (including ours) have two advantages: (1) VMR-based methods use non-interventional one-time data. Time series data require measuring the same cell multiple times without killing it, and interventional data require some techniques to interfere with gene expression, such as gene knockdown. Therefore, non-interventional one-time data used in VMR-based methods are much easier and cheaper to obtain. (2) VMR-based methods do not estimate parameters, and only calculate the mean and variance of the expression level. Some other methods need to estimate many parameters or approximate some complicated quantities, meaning that they need large data size and high data accuracy. Therefore, our method is easy to calculate, and need lower data accuracy and smaller data size.
Compared to other VMR-based methods, our method has few restrictions on the model, making them applicable to various scenarios with different dynamics. Besides, our deriva-tions do not use any approximations.
In sum, compared to other VMR-based methods, our method is universal. Compared to other non-VMR-based methods, our method is simple, and has lower requirements on data quality.
Compared to other non-VMR-based methods, our method has some disadvantages: (1) The GRN structure needs to be known. (2) Our method does not work for certain genes, depending on regulatory relationships. Proposition 1 only works for a gene that is not regulated by other genes, and we require its expression to be autonomous; Proposition 2 only works for a gene that is not in a feedback loop. (3) Proposition 2 requires the per molecule degradation rate to be a constant, and it cannot provide information about autoregulation if VMR ≥ 1. (4) Our method only works for cells at equilibrium. Thus time series data that contain time-specific information cannot be utilized other than treated as one-time data. With just the stationary distribution, sometimes it is impossible to build the causal relationship (including autoregulation) [85]. Thus with this data type, some disadvantages are inevitable. Such impossibility results might be generalized to other data types or even other fields [77].

Scenario of a single isolated gene 4.1 Setup
We first consider the expression level (e.g., mRNA count) of one gene V in a single cell. At the single-cell level, gene expression is essentially stochastic, and we do not further consider differential equation approaches [80] dynamical system approaches with deterministic [82] or stochastic [95] operators. We use a random variable X to represent the mRNA count of V . We assume V is not in a feedback loop. We also assume all environmental factors and other genes that can affect X are kept at constant levels, so that we can focus on V alone. This can be achieved if no other genes point to gene V in the GRN, such as PIP3 in Fig. 1. Then we assume that the expression of V is autonomous, thus X satisfies a time-homogeneous Markov chain defined on Z * .
Assume that the mRNA synthesis rate at X(t) = n − 1, namely the transition rate from X = n−1 to X = n, is f n ≥ 0. Assume that with n mRNA molecules, the degradation rate for each mRNA molecule is g n > 0. Then the overall degradation rate at X(t) = n, namely the transition rate from X = n to X = n − 1, is g n n. The associated master equation is (1) When f n , g n take specific forms, this master equation also corresponds to a branching process, so that related techniques can be applied [40]. Define the relative growth rate h n = f n /g n . We assume that as time tends to infinity, the process reaches equilibrium, where (1) the stationary probability distribution P n = lim t→∞ P[X(t) = n] exists, and P n = P n−1 h n /n; (2) the mean lim t→∞ E[X(t)] and the variance lim t→∞ σ 2 [X(t)] are finite. Such requirements can be satisfied under simple assumptions, such as assuming h n has a finite upper bound [53,83].
If h n > h n−1 for some n, then there exists positive autoregulation. If h n < h n−1 for some n, then there exists negative autoregulation. If there is no autoregulation, then h n is a constant h, and the stationary distribution is Poissonian with parameter h. In this setting, positive autoregulation and negative autoregulation might coexist, meaning that h n+1 < h n for some n and h n ′ +1 > h n ′ for some n ′ .

Theoretical results
With single-cell non-interventional one-time gene expression data for one gene, we have the stationary distribution of the Markov chain X. We can infer the existence of autoregulation with the VMR of X, defined as VMR(X) = σ 2 (X)/E(X). The idea is that if we let f n increase/decrease with n, and control g n to make E(X) invariant, then the variance σ 2 (X) increases/decreases [76, Section 2.5.1]. We shall prove that VMR > 1 implies the occurrence of positive autoregulation, and VMR < 1 implies the occurrence of negative autoregulation. Notice that VMR > 1 does not exclude the possibility that negative autoregulation exists for some expression level. This also applies to VMR < 1 and positive autoregulation.
We can illustrate this result with a linear model: Lemma 1. Consider a Markov chain X(t) that follows Eq. 1 with general transition coefficients f n , g n . Here X(t) models the mRNA/protein count of one gene whose expression is autonomous. Calculate VMR(X) at stationarity. (1) Assume h n+1 ≥ h n for all n. We have VMR(X) ≥ 1; moreover, VMR(X) = 1 if and only if h n+1 = h n for all n.
(2) Assume h n+1 ≤ h n for all n. We have VMR(X) ≤ 1; moreover, VMR(X) = 1 if and only if h n+1 = h n for all n.
We can take negation of Lemma 1 to obtain the following proposition.
Proposition 1. In the setting of Lemma 1, (1) If VMR(X) > 1, then there exists at least one value of n for which h n+1 > h n ; thus this gene has positive autoregulation. (2) If VMR(X) < 1, then there exists at least one value of n for which h n+1 < h n ; thus this gene has negative autoregulation. (3) If VMR(X) = 1, then either (A) h n+1 = h n for all n, meaning that this gene has no autoregulation; or (B) h n+1 < h n for one n and h n ′ +1 > h n ′ for another n ′ , meaning that this gene has both positive and negative autoregulation (at different expression levels).

Remark 1.
Results similar to Proposition 1 have been proven by Jia et al. in another model of expression for a single gene [38]. However, they require that g i = g j for any i, j. Proposition 1 can handle arbitrary g i , thus being novel.
Proof of Lemma 1. Define λ = − log P 0 , so that P 0 = exp(−λ). Define d n = n i=1 h i > 0 and stipulate that d 0 = 1. We can see that Also, Then Besides, (1) Assume h n+1 ≥ h n for all n. Then Here the first inequality is from the Cauchy-Schwarz inequality, and the second inequality is from d n d n+2 ≥ d 2 n+1 for all n. Then VMR( Equality holds if and only if d n /d n+2 = d n+1 /d n+3 for all n (the first inequality of Eq. 2) and d n d n+2 = d 2 n+1 for all n (the second inequality of Eq. 2). The equality condition is equivalent to h n+1 = h n for all n.
(2) Assume h n+1 ≤ h n for all n. Then d n+2 /d n+1 ≤ d n+1 /d n , and d n ≤ h n 1 for all n. Define Since 0 < d n ≤ h n 1 , this series converges for all t ∈ C, so that H(t) is a well-defined analytical function on C, and In the following, we only consider H(t), H ′ (t), H ′′ (t) as real functions for t ∈ R.
To prove VMR(X) ≤ 1, we just need to prove E( is a fixed interval in R with 0 < a < 1 and 1 < b < ∞. Thus t = 1 is an interior point of I. Since H(t), H ′ (t), H ′′ (t) have positive lower bounds on I, the following statements are obviously equivalent: To prove (i), we just need to prove (iv).
Consider any t 1 , t 2 ∈ I with t 1 ≤ t 2 and any p, q ∈ N with p ≥ q. Since d p+1 /d p ≤ d q+1 /d q , and t p−q 2 for all t ∈ I, and VMR(X) ≤ 1. About the condition for the equality to hold, assume h n ′ +1 < h n ′ for a given n ′ . Then for all t 1 , t 2 ∈ I with t 1 ≤ t 2 and a constant C that does not depend on t 1 , t 2 . Therefore, Since H(t) has a finite positive upper bound A and a positive lower bound B on I, we have We have proved in (1) that if h n+1 = h n for all n, then VMR(X) = 1. Thus when h n+1 ≤ h n for all n, VMR(X) = 1 if and only if h n+1 = h n for all n.
In sum, for the Markov chain model of one gene (by assuming the expression to be autonomous), when we have the stationary distribution from single-cell non-interventional one-time gene expression data, we can calculate the VMR of X. VMR(X) > 1 means the existence of positive autoregulation (while negative autoregulation might still be possible at different expression levels), and VMR(X) < 1 means the existence of negative autoregulation (while positive autoregulation might still be possible at different expression levels). VMR(X) = 1 means either (1) no autoregulation exists; or (2) both positive autoregulation and negative autoregulation exist (at different expression levels). In reality, many genes are non-autonomous, and transcriptional/translational bursting can make the VMR to be larger than 100 [54]. Since Proposition 1 does not apply to non-autonomous cases, such genes might not have autoregulations.

Setup
We consider m genes V 1 , . . . , V m for a single cell. Denote their expression levels by random variables X 1 , . . . , X m . The change of X i can depend on X j (mutual regulation) and X i itself (autoregulation). Since these genes regulate each other, and their expression levels are not fixed, we cannot consider them separately. If the expression of gene V k is nonautonomous, we also need to add its interior factors (gene state and/or mRNA count) into X 1 , . . . , X m .
We can use a continuous-time one-step Markov chain on (Z * ) m to describe the dynamics. Each state of this Markov chain, (X 1 = n 1 , . . . , X i = n i , . . . , X m = n m ), can be abbreviated as n = (n 1 , . . . , n i , . . . , n m ). For gene V i , the transition rate of n i − 1 → n i is f i (n), and the transition rate of n i → n i − 1 is g i (n)n i . Transitions with more than one step are not allowed. The master equation of this process is dP(n) dt = i P(n 1 , . . . , n i + 1, . . . , n m )g i (n 1 , . . . , n i + 1, . . . , n m )(n i + 1) Define nī = (n 1 , . . . , n i−1 , n i+1 , . . . , n m ). Define h i (n) = f i (n)/g i (n) to be the relative growth rate of gene V i . Autoregulation means for some fixed nī, h i (n) is (locally) increasing/decreasing with n i , thus f i (n) increases/decreases and/or g i (n) decreases/increases with n i . For the non-autonomous scenario, another possibility for autoregulation is that V i can affect its interior factors (gene state and/or mRNA count).

Theoretical results
With expression data for multiple genes, there are various methods to infer the regulatory relationships between different genes, so that the GRN can be reconstructed [86]. In the GRN, if there is a directed path from gene V i to gene V j , meaning that V i can directly or indirectly regulate V j , then V i is an ancestor of V j , and V j is a descendant of V i . Fix a gene V k in a GRN. We consider a simple case that V k is not contained in any directed cycle (feedback loop), which means no gene is both an ancestor and a descendant of V k , such as PIP2 in Fig. 1. This means V k itself is a strongly connected component of the GRN. This condition is automatically satisfied if the GRN has no directed cycle. If the expression of V k is non-autonomous, we need to add the interior factors (gene state and/or mRNA count) of V k into V 1 , . . . , V m , and it is acceptable that V k regulates its interior factors. In this case, if the one-step model holds, we can prove that if V k does not regulate itself, meaning that h k (n) is a constant for fixed nk and different n k , and X k does not affect its interior factors (if non-autonomous), then VMR(X k ) ≥ 1. The reason is that VMR < 1 requires either a feedback loop or autoregulation. Certainly, VMR < 1 might also mean that the one-step model fails. One intuition is to assume the transitions of Vk are extremely slow, so that V k is approximately the average of many Poisson variables. It is easy to verify that the average of Poisson variables has VMR ≥ 1. We need to assume that the per molecule degradation rate g k (·) for V k is not affected by V 1 , . . . , V m , which is not always true in reality [42]. With this result, when VMR < 1, there might be autoregulation.

Proposition 2.
Consider the one-step Markov chain model for multiple genes, described by Eq. 3. Assume the GRN has no directed cycle, or at least there is no directed cycle that contains gene V k . Assume g k (·) is a constant for all n. If V k has no autoregulation, meaning that h k (·) and f k (·) do not depend on n k , and V k does not regulate its interior factors (gene state and/or mRNA count), then V k has VMR ≥ 1. Therefore, V k has VMR < 1 means V k has autoregulation, or the one-step model fails. [29,93], and they state Proposition 2 in an unpublished work. Proposition 2 also appears in a preprint by Mahajan et al. [49], but the proof is based on a linear noise approximation, which requires that f k (·) is linear with nk. We propose a rigorous proof independently.

Paulsson et al. study a similar problem
Proof. Denote the expression level of V k by W . Assume the ancestors of V k are V 1 , . . . , V l . For simplicity, denote the expression levels of V 1 , . . . , V l by a (high-dimensional) random variable Y . Assume V k has no autoregulation. Since V k does not regulate V 1 , . . . , V l , W does not affect Y . Denote the transition rate from Y = i to Y = j by q ij ≥ 0. Stipulate that q ii = − j =i q ij . When Y = i, the transition rate from W = n to W = n + 1 is F i (does not depend on n), and the transition rate from W = n to W = n − 1 is G.
The master equation of this process is Assume there is a unique stationary probability distribution P n,i = lim t→∞ P[W (t) = n, Y (t) = i]. This can be guaranteed by assuming the process to be irreducible. Then we have Define P i = n P n,i . Sum over n for Eq. 4 to obtain meaning that P i is the stationary probability distribution of Y .
Define W i to be W conditioned on Y = i at stationarity. Then P(W i = n) = P(W = n | Y = i) = P n,i /P i , and E(W i ) = n nP n,i /P i . Multiply Eq. 4 by n and sum over n to obtain Here and in the following, we repeatedly apply the tricks of splitting n and shifting the index of summation. For example, Sum over i for Eq. 6 to obtain Multiply Eq. 4 by n 2 and sum over n to obtain Sum over i for Eq. 8 to obtain Multiply Eq. 6 by E(W i ) and sum over i to obtain Then we have Here the first equality is from Eq. 10, the third equality is from Eq. 5, and other equalities are equivalent transformations. Now we have where the first equality is by definition, the second equality is from Eqs. 7,9, the first inequality is from Eq. 11, the third equality is from i P i = 1, and the second inequality is the Cauchy-Schwarz inequality.

Remark 2.
In gene expression, the total noise (σ 2 (X)/(EX) 2 ) can be decomposed into intrinsic (cellular) noise and extrinsic (environmental) noise [3,71,27,47,75]. Inspired by that, we can decompose the VMR into intrinsic and extrinsic components. Denote intrinsic and extrinsic stochastic factors as I, E, and the expression level X is a deterministic function of these factors: X = X(I, E). Then where E I|E is the expectation conditioned on E. This decomposition might lead to further understanding of Proposition 2.
We hypothesize that the requirement for g k (·) in Proposition 2 can be dropped: Assume V k is not contained in a directed cycle in the GRN, and V k does not regulate its interior factors (gene state and/or mRNA count). If V k has no autoregulation, meaning that h k (·) does not depend on n k (but might depend on nk), then V k has VMR ≥ 1.
The main obstacle of proving this conjecture is that the second equality in Eq. 12 does not hold. The reason is that G i cannot be extracted from the summation, and we cannot link i P i E(W 2 i ) and i G i P i E(W 2 i ). If the GRN has directed cycles, there is a result by Paulsson et al. [29,93], which is proved under first-order approximations of covariances. The general case (when the approximations do not apply) has been numerically verified but not proved yet: does not depend on n, and f i (·) does not depend on n i (no autoregulation). Then for at least one gene V j , we have VMR ≥ 1 [29,93].
Due to the existence of directed cycles, one gene can affect itself through other genes, and we cannot study them separately.
Assume Conjecture 2 is correct. For m genes, if we find that VMR for each gene is less than 1, then we can infer that autoregulation exists, although we do not know which gene has autoregulation. Another possibility is that the one-step model fails.

Applying theoretical results to experimental data
We summarize our theoretical results into Algorithm 1. Proposition 1 applies to a gene that has no ancestor in the GRN. However, it requires the corresponding gene has autonomous expression (or the transition rates of gene states are high enough, so that the non-autonomous process is close to an autonomous process), which is difficult to validate and often does not hold in reality. Thus the inference result by Proposition 1 for VMR > 1 (positive autoregulation) is not very reliable. When VMR < 1 and Proposition 1 could apply, we should instead apply Proposition 2 to determine the existence of autoregulation, since Proposition 2 does not require the expression to be autonomous, thus being much more reliable, although it may fail if the one-step model does not hold. Proposition 2 applies when the gene is not in a feedback loop and has VMR < 1. Notice that our result cannot determine that a gene has no autoregulation.
For a given gene without autoregulation, its expression level satisfies a Poisson distribution, and VMR is 1. If we have n samples of its expression level, then the sample VMR (sample variance divided by sample mean) asymptotically satisfies a Gamma distribution Γ[(n−1)/2, 2/(n−1)], and we can determine the confidence interval of sample VMR [20]. If the sample VMR is out of this confidence interval, then we know that VMR is significantly different from 1, and Propositions 1,2 might apply.
We apply our method to four groups of single-cell non-interventional one-time gene expression data from experiments, where the corresponding GRNs are known. Notice that we need to convert indirect measurements into protein/mRNA count. See Table 1 for our inference results and theoretical/experimental evidence that partially validates our results. See Appendix C for details. There are 186 genes in these four data sets, and we can only determine that 12 genes have autoregulation (7 genes determined by Proposition 1, and 5 genes determined by Proposition 2). Not every VMR is less than 1, so that Conjecture 2 does not apply. For the other 174 genes, (1) some of them are not contained in the known GRN, and we cannot determine if they are in directed cycles; (2) some of them are in directed cycles; (3) some of them have ancestors, and we cannot reject the hypothesis that VMR ≥ 1; (4) some of them have no ancestors, and we cannot reject the hypothesis that VMR = 1. Therefore, Proposition 1 and Proposition 2 do not apply, and we do not know whether they have autoregulation.
In some cases, we have experimental evidence that some genes have autoregulation, so that we can partially validate our inference results. Nevertheless, as discussed in the Introduction, there is no gold standard to evaluate our inference results. Besides, Proposition 2 requires that the one-step model holds, which we cannot verify.
In the data set by Guo et al. [26], Sanchez-Castillo et al. [59] inferred that 17 of 39 genes have autoregulation, and 22 genes do not have autoregulation. We infer that 5 genes have autoregulation, and 34 genes cannot be determined. Here 3 genes are shared by both inference results to have autoregulation. Consider a random classifier that randomly picks 5 genes and claims they have autoregulation. Using Sanchez-Castillo et al. as the standard, this random classifier has probability 62.55% to be worse than our result, and 10.17% to be better than our result. Thus our inference result is better than a random classifier, but the advantage is not substantial.  Table 1: The autoregulation inference results by our method on four data sets. Source column is the paper that contains this data set. Proposition 1 column is the genes that can be only inferred by Proposition 1 to have autoregulation. Proposition 2 column is the genes that can be inferred by Proposition 2 to have autoregulation. Theory column is the genes inferred by both our method and other theoretical works to have autoregulation. Experiment column is the genes inferred by both our method and other experimental works to have autoregulation. Bold font means the inferred gene with autoregulation is validated by other results. Details can be found in Appendix C.

Conclusions
For a single gene that is not affected by other genes, or a group of genes that form a connected GRN, we develop rigorous theoretical results (without applying approximations) to determine the existence of autoregulation. These results generalize known relationships between autoregulation and VMR by dropping restrictions on parameters. Our results only depend on VMR, which is easy to compute and more robust than other complicated statistics. We also apply our method to experimental data and detect some genes that might have autoregulation. Our method requires independent and identically distributed samples from the exact stationary distribution of a fully observed Markov chain, plus a known GRN. Proposition 1 requires that the expression is autonomous. Proposition 2 requires that the Markov chain model is one-step, the GRN has no directed cycle, and degradation is not regulated. If our inference fails, then some requirements are not met: (1) cells might affect each other, making the samples dependent; (2) cells are heterogeneous; (3) the measurements have extra errors; (4) the cells are not at stationarity; (5) there exist unobserved variables that affect gene expression; (6) the GRN is inferred by a theoretical method, which can be interfered by the existence of autoregulation; (7) the expression is non-autonomous; (8) the Markov chain is multi-step; (9) the GRN has unknown directed cycles; (10) the degradation rate is regulated by other genes. Such situations, especially the unobserved variables, are unavoidable. Therefore, current data might not satisfy these requirements, and our inference results should be interpreted as informative findings, not ground truths.
There are some known methods that overcome the above obstacles, and there are also some possible solutions that might appear in the future. (1) The dependency can be solved by better measurements for isolated cells that do not affect each other. In fact, the relationship between autoregulation and cell-cell interaction has been studied [46].
(2) About cell heterogeneity, we prove a result in Appendix D that if several cell types have VMR ≥ 1, then for a mixed population of such cell types, we still have VMR ≥ 1. Therefore, cell heterogeneity does not fail Proposition 2, since VMR < 1 for the mixture of several cell types means VMR < 1 for at least one cell type. (3) With the development of experimental technologies, we expect that the measurement error can decrease. (4) Some works study autoregulation in non-stationary situations [10,69,66,36]. (5) Since hidden variables hurt any mechanism-based models, we can develop methods (especially with machine learning tools) that determine autoregulation based on similarities between gene expression profiles [87,94,78,88,89]. (6) Some GRN inference methods can also determine the existence of autoregulation [59]. (7) Many methods (including our Proposition 2) work in non-autonomous situations. (8) Some works study multi-step models [7,43,74]. (9) We expect the appearance of more advanced GRN inference methods. (10) If probabilists can prove Conjecture 1, then the restriction on degradation rate can be lifted.
In fact, other theoretical works that determine gene autoregulation, or general gene

Declarations
• Funding: This research was partially supported by NIH grant R01HL146552 (Y.W.).
• Competing interests: The authors declare no conflict of interest.
• Data and code availability: All code files are available in https://github.com/YueWangMathbio/Autoregulation.

A.1 Test Proposition 1 without autoregulation
Simulation 1: Consider the Markov chain in Example 1 with k = 1, b = 0, c = 1. The stationary distribution is Poissonian with parameter 1. This process has no autoregulation with the true VMR = 1. For sample sizes n = 100, n = 1000, n = 10000, we repeat the experiment for 10000 times and calculate the rate that the sample VMR falls in the 95% confidence interval. See Table 2 for results. Proposition 1 has about 95% probability to produce the correct result that VMR = 1 (no autoregulation), since the confidence interval is 95%.
The stationary distribution is geometric with parameter 0.5. This process has positive autoregulation with the true VMR = 2. For sample sizes n = 100, n = 1000, n = 10000, we repeat the experiment for 10000 times and calculate the rate that the sample VMR falls in the 95% confidence interval. See Table 3 for results. When n is not too small, Proposition 1 always produces the correct result that VMR > 1 (positive autoregulation). For sample sizes n = 100, n = 1000, n = 10000, we repeat the experiment for 10000 times and calculate the rate that the sample VMR falls in the 95% confidence interval. See Table 4 for results. When n increases, we are very likely to obtain the correct VMR > 1, but Proposition 2 cannot determine whether autoregulation exists. 56.6% 43.4% 0% n = 1000 100% 0% 0% n = 10000 100% 0% 0% Example 1 with k = 10, b = −1, c = 1. Restricted on G = 1, M is the same as Example 1 with k = 10, b = −1, c = 2. This process has autoregulation with the true VMR = 0.733. For sample sizes n = 100, n = 1000, n = 10000, we repeat the experiment for 10000 times and calculate the rate that the sample VMR falls in the 95% confidence interval. See Table 5 for results. When n is not too small, Proposition 2 always produces the correct result that VMR < 1 (autoregulation).
After normalization, we find that c = 2e/(4e − 1). Thus . Besides, Then for the first gene X, we have Thus Due to symmetry, the other gene also has VMR(Y ) ≈ 0.55.

C Details of applications on experimental data
In experiments, the expression levels of genes are not directly measured as mRNA or protein counts. Rather, they are measured as cycle threshold (Ct) values or fluorescence intensity values. Such indirect measurements need to be converted. Related details can be found in other papers [38]. Guo et al. [26] measured the expression (mRNA) levels of 48 genes for mouse embryo cells at different developmental stages. We consider three groups (16-cell stage, 32-cell stage, 64-cell stage) that have more than 50 samples. Sanchez-Castillo et al. [59] used such data to infer the GRN structure, including autoregulation, but the inferred GRN only contains 39 genes. We cannot guarantee that the other 9 genes have no ancestors in the true GRN (to apply Proposition 1) or these genes are not contained in directed cycles (to apply Proposition 2). Thus we ignore those 9 genes not in this GRN. In the inferred GRN, genes BMP4, CREB312, and TCFAP2C are not contained in directed cycles. In the 16- Thus we can apply Proposition 1 to infer that BIM, CCND1, ECT2, PFKP might have positive autoregulation. Nevertheless, it is more likely that the expressions of these four genes are non-autonomous, and there is no autoregulation. There is experimental evidence that ECT2 has autoregulation [28], which partially validates our inference results. No other gene fits the requirement of Proposition 2.
Moignard et al. [50] measured the expression (mRNA) levels of 46 genes for mouse embryo cells. Chan et al. [13] inferred the GRN structure (autoregulation not included). Gene EIF2B1 has 3934 effective samples, and VMR is 0.66, outside of the 95% confidence interval [0.9563, 1.0447]. Gene EIF2B1 has 12 effective samples, and VMR is 0.24, outside of the 95% confidence interval [0.3469, 1.9927]. We can apply Proposition 2 to infer that EIF2B1 and HOXD8 might have autoregulation. No other gene fits the requirement of Proposition 1.
Sachs et al. [58] measured the expression (protein) levels of 11 genes in the RAF signaling pathway for human T cells. The measurements were repeated for 14 groups of cells under different interventions. Werhli et al. [91] inferred the GRN structure (autoregulation not included). In the inferred GRN (Fig. 1), PIP3 gene has no ancestor, and its VMRs in all 14 groups are larger than 5, while the 95% confidence intervals for all 14 groups are contained in [0.8, 1.2]. Therefore, we can apply Proposition 1 and infer that PIP3 might have positive autoregulation. Nevertheless, it is more likely that the expression of PIP3 is non-autonomous, and there is no autoregulation. No other gene fits the requirement of Proposition 2.

D Heterogeneity and VMR
Proposition 3. Consider n independent random variables X 1 , . . . , X n and probabilities p 1 , . . . , p n with p i = 1. Consider an independent random variable R that equals i with probability p i . Construct a random variable Z that equals X i when R = i. If each X i has VMR ≥ 1, then Z has VMR ≥ 1.
Proof. We only need to prove this for n = 2. The case for general n can be proved by mixing two variables iteratively.
Consider random variables X, Y and construct Z that equals X or Y with probability p or 1 − p. Since VMR(X) ≥ 1, VMR(Y ) ≥ 1, we have E(X 2 ) − [E(X)] 2 ≥ E(X) and